simulation results/results_pois_data/ratio of lambdas sim res.R

In BriceonWiley/IntegratedLikelihood.R:

library(tidyverse); theme_set(theme_minimal(15))
setwd("/Users/BriceonWiley/Documents/Coding/RStudio/Research/two mean missclassification/sim_results")

produce_results <- function(data, type) {
  type <- tolower(type)
  if (type == "coverage") map_dfr(data, ~ coverage_phi(.x, TRUE))
  else if (type == "width") map_dfr(data, ~ width_phi(.x, TRUE))
}
convert_coverage <- function(data, adj) {
  if (adj) {
    data %>%
      mutate(
        W = case_when(
          (W >= 0.98) ~ "c > 98%",
          (W >= 0.96 & W < 0.98) ~ "96% < c < 98%",
          (W >= 0.94 & W < 0.96) ~ "94% < c < 96%",
          (W >= 0.92 & W < 0.94) ~ "92% < c < 94%",
          (W >= 0.90 & W < 0.92) ~ "90% < c < 92%",
          (W < 0.90) ~ "c < 90%"
        ),
        S = case_when(
          (S >= 0.98) ~ "c > 98%",
          (S >= 0.96 & S < 0.98) ~ "96% < c < 98%",
          (S >= 0.94 & S < 0.96) ~ "94% < c < 96%",
          (S >= 0.92 & S < 0.94) ~ "92% < c < 94%",
          (S >= 0.90 & S < 0.92) ~ "90% < c < 92%",
          (S < 0.90) ~ "c < 90%"
        ),
        ILR = case_when(
          (ILR >= 0.98) ~ "c > 98%",
          (ILR >= 0.96 & ILR < 0.98) ~ "96% < c < 98%",
          (ILR >= 0.94 & ILR < 0.96) ~ "94% < c < 96%",
          (ILR >= 0.92 & ILR < 0.94) ~ "92% < c < 94%",
          (ILR >= 0.90 & ILR < 0.92) ~ "90% < c < 92%",
          (ILR < 0.90) ~ "c < 90%"
        )
      )
  } else {
    data %>%
      mutate(
        W = W - 0.95,
        S = S - 0.95,
        ILR = ILR - 0.95
      )
  }
}
coverage_plots <- function(data) {
  data %>%
    mutate(
      W = case_when(
        (W >= 0.98) ~ "c > 98%",
        (W >= 0.96 & W < 0.98) ~ "96% < c < 98%",
        (W >= 0.94 & W < 0.96) ~ "94% < c < 96%",
        (W >= 0.92 & W < 0.94) ~ "92% < c < 94%",
        (W >= 0.90 & W < 0.92) ~ "90% < c < 92%",
        (W < 0.90) ~ "c < 90%"
      ),
      S = case_when(
        (S >= 0.98) ~ "c > 98%",
        (S >= 0.96 & S < 0.98) ~ "96% < c < 98%",
        (S >= 0.94 & S < 0.96) ~ "94% < c < 96%",
        (S >= 0.92 & S < 0.94) ~ "92% < c < 94%",
        (S >= 0.90 & S < 0.92) ~ "90% < c < 92%",
        (S < 0.90) ~ "c < 90%"
      ),
      ILR = case_when(
        (ILR >= 0.98) ~ "c > 98%",
        (ILR >= 0.96 & ILR < 0.98) ~ "96% < c < 98%",
        (ILR >= 0.94 & ILR < 0.96) ~ "94% < c < 96%",
        (ILR >= 0.92 & ILR < 0.94) ~ "92% < c < 94%",
        (ILR >= 0.90 & ILR < 0.92) ~ "90% < c < 92%",
        (ILR < 0.90) ~ "c < 90%"
      )
    ) %>%
    select(-N0, -N, -theta2) %>%
    gather("Interval", "Coverage", -phi, -theta1) %>%
    mutate(
      Coverage = fct_relevel(
        Coverage, c(
          "c > 98%", "96% < c < 98%", "94% < c < 96%",
          "92% < c < 94%", "90% < c < 92%", "c < 90%"
        )
      ),
      Interval = fct_relevel(Interval, c("W", "ILR", "S"))
    ) %>%
    ggplot(aes(theta1, phi)) +
    geom_tile(aes(fill = Coverage)) +
    facet_wrap(~ Interval) +
    scale_fill_manual(
      values = c(
        "c > 98%" = "#053061",
        "96% < c < 98%" = "#4393c3",
        "94% < c < 96%" = "#f7f7f7",
        "92% < c < 94%" = "#f4a582",
        "90% < c < 92%" = "#b2182b",
        "c < 90%" = "#67001f"
      )
    ) +
    scale_x_continuous(
      name = expression(theta[1]==theta[2]),
      breaks = c(0.05, 0.50, 0.95)
    ) +
    labs(y = expression(phi), fill = "Coverage\nPerformance") +
    theme(
      axis.title.y = element_text(angle = 0, vjust = 0.5, size = 18),
      axis.title.x = element_text(size = 18)
    )
}
width_plots <- function(data) {
  data %>%
    select(-N0, -N, -theta2) %>%
    gather("Interval", "Width", -phi, -theta1) %>%
    mutate(
      Interval = fct_relevel(Interval, c("W", "ILR", "S"))
    ) %>%
    ggplot(aes(theta1, phi)) +
    geom_tile(aes(fill = Width)) +
    facet_wrap(~ Interval) +
    scale_fill_gradient(
      low = "white",
      high = "black"
    ) +
    scale_x_continuous(
      name = expression(theta[1]==theta[2]),
      breaks = c(0.05, 0.50, 0.95)
    ) +
    labs(y = expression(phi), fill = "Average\nWidth") +
    theme(
      axis.title.y = element_text(angle = 0, vjust = 0.5, size = 18),
      axis.title.x = element_text(size = 18)
    )
}


############################################### Code for accessing results ####
# sim_regular <- list(
#   "N0_1_N_2" = phi_N0_1_N_2,
#   "N0_1_N_3" = phi_N0_1_N_3,
#   "N0_2_N_4" = phi_N0_2_N_4,
#   "N0_2_N_6" = phi_N0_2_N_6,
#   "N0_3_N_6" = phi_N0_3_N_6,
#   "N0_3_N_9" = phi_N0_3_N_9,
#   "N0_4_N_8" = phi_N0_4_N_8
# )
# res_cov <- map(sim_regular, ~ produce_results(.x, "coverage"))
# res_wid <- map(sim_regular, ~ produce_results(.x, "width"))
# saveRDS(res_cov, "res_cov.rds")
# saveRDS(res_wid, "res_wid.rds")
res_cov <- readRDS("res_cov.rds")
res_wid <- readRDS("res_wid.rds")
###############################################################################

figs_cov <- map(res_cov, coverage_plots)
figs_wid <- map(res_wid, width_plots)

map2(
  figs_cov, names(figs_cov),
  ~ ggsave(glue("c04_cov_{.y}.pdf"), .x)
)
map2(
  figs_wid, names(figs_cov),
  ~ ggsave(glue("c04_wid_{.y}.pdf"), .x)
)

################################################## For combo width figures ####
wid_2_6_3_6 <-
  bind_rows(res_wid$N0_2_N_6, res_wid$N0_3_N_6, .id = "group") %>%
  mutate(
    group = case_when(
      (group == 1) ~ "(2,6)",
      (group == 2) ~ "(3,6)"
    ) %>%
      as_factor()
  ) %>%
  select(-N0, -N, -theta2) %>%
  gather("Interval", "Width", -phi, -theta1, -group) %>%
  mutate(
    Interval = fct_relevel(Interval, c("W", "ILR", "S"))
  ) %>%
  ggplot(aes(theta1, phi)) +
  geom_tile(aes(fill = Width)) +
  facet_grid(rows = vars(group), cols = vars(Interval)) +
  scale_fill_gradient(
    low = "white",
    high = "black"
  ) +
  scale_x_continuous(
    name = expression(theta[1]),
    breaks = c(0.05, 0.50, 0.95)
  ) +
  labs(y = expression(phi), fill = "Average\nWidth") +
  theme(
    axis.title.y = element_text(angle = 0, vjust = 0.5),
    strip.text.y = element_text(angle=0)
  )

wid_3_9_4_8 <-
  bind_rows(res_wid$N0_3_N_9, res_wid$N0_4_N_8, .id = "group") %>%
  mutate(
    group = case_when(
      (group == 1) ~ "(3,9)",
      (group == 2) ~ "(4,8)"
    ) %>%
      as_factor()
  ) %>%
  select(-N0, -N, -theta2) %>%
  gather("Interval", "Width", -phi, -theta1, -group) %>%
  mutate(
    Interval = fct_relevel(Interval, c("W", "ILR", "S"))
  ) %>%
  ggplot(aes(theta1, phi)) +
  geom_tile(aes(fill = Width)) +
  facet_grid(rows = vars(group), cols = vars(Interval)) +
  scale_fill_gradient(
    low = "white",
    high = "black"
  ) +
  scale_x_continuous(
    name = expression(theta[1]),
    breaks = c(0.05, 0.50, 0.95)
  ) +
  labs(y = expression(phi), fill = "Average\nWidth") +
  theme(
    axis.title.y = element_text(angle = 0, vjust = 0.5),
    strip.text.y = element_text(angle=0)
  )

wid_2_6_3_6_3_9_4_8 <-
  bind_rows(
    res_wid$N0_2_N_6,
    res_wid$N0_3_N_6,
    res_wid$N0_3_N_9,
    res_wid$N0_4_N_8,
    .id = "group"
  ) %>%
  mutate(
    group = case_when(
      (group == 1) ~ "(2,6)",
      (group == 2) ~ "(3,6)",
      (group == 3) ~ "(3,9)",
      (group == 4) ~ "(4,8)"
    ) %>%
      as_factor()
  ) %>%
  select(-N0, -N, -theta2) %>%
  gather("Interval", "Width", -phi, -theta1, -group) %>%
  mutate(
    Interval = fct_relevel(Interval, c("W", "ILR", "S"))
  ) %>%
  ggplot(aes(theta1, phi)) +
  geom_tile(aes(fill = Width)) +
  facet_grid(rows = vars(group), cols = vars(Interval)) +
  scale_fill_gradient(
    low = "white",
    high = "black"
  ) +
  scale_x_continuous(
    name = expression(theta[1]==theta[2]),
    breaks = c(0.05, 0.50, 0.95)
  ) +
  labs(y = expression(phi), fill = "Average\nWidth") +
  theme(
    axis.title.y = element_text(angle = 0, vjust = 0.5, size = 18),
    axis.title.x = element_text(size = 18)
  )
############################################### For combo coverage figures ####
cov_2_6_3_6 <-
  bind_rows(res_cov$N0_2_N_6, res_cov$N0_3_N_6, .id = "group") %>%
  mutate(
    group = case_when(
      (group == 1) ~ "(2,6)",
      (group == 2) ~ "(3,6)"
    ) %>%
      as_factor()
  ) %>%
  mutate(
    W = case_when(
      (W >= 0.98) ~ "c > 98%",
      (W >= 0.96 & W < 0.98) ~ "96% < c < 98%",
      (W >= 0.94 & W < 0.96) ~ "94% < c < 96%",
      (W >= 0.92 & W < 0.94) ~ "92% < c < 94%",
      (W >= 0.90 & W < 0.92) ~ "90% < c < 92%",
      (W < 0.90) ~ "c < 90%"
    ),
    S = case_when(
      (S >= 0.98) ~ "c > 98%",
      (S >= 0.96 & S < 0.98) ~ "96% < c < 98%",
      (S >= 0.94 & S < 0.96) ~ "94% < c < 96%",
      (S >= 0.92 & S < 0.94) ~ "92% < c < 94%",
      (S >= 0.90 & S < 0.92) ~ "90% < c < 92%",
      (S < 0.90) ~ "c < 90%"
    ),
    ILR = case_when(
      (ILR >= 0.98) ~ "c > 98%",
      (ILR >= 0.96 & ILR < 0.98) ~ "96% < c < 98%",
      (ILR >= 0.94 & ILR < 0.96) ~ "94% < c < 96%",
      (ILR >= 0.92 & ILR < 0.94) ~ "92% < c < 94%",
      (ILR >= 0.90 & ILR < 0.92) ~ "90% < c < 92%",
      (ILR < 0.90) ~ "c < 90%"
    )
  ) %>%
  select(-N0, -N, -theta2) %>%
  gather("Interval", "Coverage", -phi, -theta1, -group) %>%
  mutate(
    Coverage = fct_relevel(
      Coverage, c(
        "c > 98%", "96% < c < 98%", "94% < c < 96%",
        "92% < c < 94%", "90% < c < 92%", "c < 90%"
      )
    ),
    Interval = fct_relevel(Interval, c("W", "ILR", "S"))
  ) %>%
  ggplot(aes(theta1, phi)) +
  geom_tile(aes(fill = Coverage)) +
  facet_grid(rows = vars(group), cols = vars(Interval)) +
  scale_fill_manual(
    values = c(
      "c > 98%" = "#053061",
      "96% < c < 98%" = "#4393c3",
      "94% < c < 96%" = "#f7f7f7",
      "92% < c < 94%" = "#f4a582",
      "90% < c < 92%" = "#b2182b",
      "c < 90%" = "#67001f"
    )
  ) +
  scale_x_continuous(
    name = expression(theta[1]),
    breaks = c(0.05, 0.50, 0.95)
  ) +
  labs(y = expression(phi), fill = "Coverage\nPerformance") +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5)) +
  theme(
    axis.title.y = element_text(angle = 0, vjust = 0.5),
    strip.text.y = element_text(angle=0)
  )

cov_3_9_4_8 <-
  bind_rows(res_cov$N0_3_N_9, res_cov$N0_4_N_8, .id = "group") %>%
  mutate(
    group = case_when(
      (group == 1) ~ "(3,9)",
      (group == 2) ~ "(4,8)"
    ) %>%
      as_factor()
  ) %>%
  mutate(
    W = case_when(
      (W >= 0.98) ~ "c > 98%",
      (W >= 0.96 & W < 0.98) ~ "96% < c < 98%",
      (W >= 0.94 & W < 0.96) ~ "94% < c < 96%",
      (W >= 0.92 & W < 0.94) ~ "92% < c < 94%",
      (W >= 0.90 & W < 0.92) ~ "90% < c < 92%",
      (W < 0.90) ~ "c < 90%"
    ),
    S = case_when(
      (S >= 0.98) ~ "c > 98%",
      (S >= 0.96 & S < 0.98) ~ "96% < c < 98%",
      (S >= 0.94 & S < 0.96) ~ "94% < c < 96%",
      (S >= 0.92 & S < 0.94) ~ "92% < c < 94%",
      (S >= 0.90 & S < 0.92) ~ "90% < c < 92%",
      (S < 0.90) ~ "c < 90%"
    ),
    ILR = case_when(
      (ILR >= 0.98) ~ "c > 98%",
      (ILR >= 0.96 & ILR < 0.98) ~ "96% < c < 98%",
      (ILR >= 0.94 & ILR < 0.96) ~ "94% < c < 96%",
      (ILR >= 0.92 & ILR < 0.94) ~ "92% < c < 94%",
      (ILR >= 0.90 & ILR < 0.92) ~ "90% < c < 92%",
      (ILR < 0.90) ~ "c < 90%"
    )
  ) %>%
  select(-N0, -N, -theta2) %>%
  gather("Interval", "Coverage", -phi, -theta1, -group) %>%
  mutate(
    Coverage = fct_relevel(
      Coverage, c(
        "c > 98%", "96% < c < 98%", "94% < c < 96%",
        "92% < c < 94%", "90% < c < 92%", "c < 90%"
      )
    ),
    Interval = fct_relevel(Interval, c("W", "ILR", "S"))
  ) %>%
  ggplot(aes(theta1, phi)) +
  geom_tile(aes(fill = Coverage)) +
  facet_grid(rows = vars(group), cols = vars(Interval)) +
  scale_fill_manual(
    values = c(
      "c > 98%" = "#053061",
      "96% < c < 98%" = "#4393c3",
      "94% < c < 96%" = "#f7f7f7",
      "92% < c < 94%" = "#f4a582",
      "90% < c < 92%" = "#b2182b",
      "c < 90%" = "#67001f"
    )
  ) +
  scale_x_continuous(
    name = expression(theta[1]),
    breaks = c(0.05, 0.50, 0.95)
  ) +
  labs(y = expression(phi), fill = "Coverage\nPerformance") +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5)) +
  theme(
    axis.title.y = element_text(angle = 0, vjust = 0.5),
    strip.text.y = element_text(angle=0)
  )

cov_2_6_3_6_3_9_4_8 <-
  bind_rows(
    res_cov$N0_2_N_6,
    res_cov$N0_3_N_6,
    res_cov$N0_3_N_9,
    res_cov$N0_4_N_8,
    .id = "group"
  ) %>%
  mutate(
    group = case_when(
      (group == 1) ~ "(2,6)",
      (group == 2) ~ "(3,6)",
      (group == 3) ~ "(3,9)",
      (group == 4) ~ "(4,8)"
    ) %>%
      as_factor()
  ) %>%
  mutate(
    W = case_when(
      (W >= 0.98) ~ "c > 98%",
      (W >= 0.96 & W < 0.98) ~ "96% < c < 98%",
      (W >= 0.94 & W < 0.96) ~ "94% < c < 96%",
      (W >= 0.92 & W < 0.94) ~ "92% < c < 94%",
      (W >= 0.90 & W < 0.92) ~ "90% < c < 92%",
      (W < 0.90) ~ "c < 90%"
    ),
    S = case_when(
      (S >= 0.98) ~ "c > 98%",
      (S >= 0.96 & S < 0.98) ~ "96% < c < 98%",
      (S >= 0.94 & S < 0.96) ~ "94% < c < 96%",
      (S >= 0.92 & S < 0.94) ~ "92% < c < 94%",
      (S >= 0.90 & S < 0.92) ~ "90% < c < 92%",
      (S < 0.90) ~ "c < 90%"
    ),
    ILR = case_when(
      (ILR >= 0.98) ~ "c > 98%",
      (ILR >= 0.96 & ILR < 0.98) ~ "96% < c < 98%",
      (ILR >= 0.94 & ILR < 0.96) ~ "94% < c < 96%",
      (ILR >= 0.92 & ILR < 0.94) ~ "92% < c < 94%",
      (ILR >= 0.90 & ILR < 0.92) ~ "90% < c < 92%",
      (ILR < 0.90) ~ "c < 90%"
    )
  ) %>%
  select(-N0, -N, -theta2) %>%
  gather("Interval", "Coverage", -phi, -theta1, -group) %>%
  mutate(
    Coverage = fct_relevel(
      Coverage, c(
        "c > 98%", "96% < c < 98%", "94% < c < 96%",
        "92% < c < 94%", "90% < c < 92%", "c < 90%"
      )
    ),
    Interval = fct_relevel(Interval, c("W", "ILR", "S"))
  ) %>%
  ggplot(aes(theta1, phi)) +
  geom_tile(aes(fill = Coverage)) +
  facet_grid(rows = vars(group), cols = vars(Interval)) +
  scale_fill_manual(
    values = c(
      "c > 98%" = "#053061",
      "96% < c < 98%" = "#4393c3",
      "94% < c < 96%" = "#f7f7f7",
      "92% < c < 94%" = "#f4a582",
      "90% < c < 92%" = "#b2182b",
      "c < 90%" = "#67001f"
    )
  ) +
  scale_x_continuous(
    name = expression(theta[1]==theta[2]),
    breaks = c(0.05, 0.50, 0.95)
  ) +
  labs(y = expression(phi), fill = "Coverage\nPerformance") +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5)) +
  theme(
    axis.title.y = element_text(angle = 0, vjust = 0.5, size = 18),
    axis.title.x = element_text(size = 18)
  )